# packages required
packages <- c("dplyr", "ggplot2", "haven", "sf", "fixest", "gt", "tidyr", "modelsummary", "patchwork", "foreign", "webshot2")

# Install packages if not already installed
installed <- packages %in% rownames(installed.packages())
if(any(!installed)) install.packages(packages[!installed])

# Load all packages
lapply(packages, library, character.only = TRUE)

# =================================================================================
# Data 
# =================================================================================
# Loading the data 

# Variable Descriptions:
# aimag                  : Aimag code (province-level administrative unit)
# YAChiglel5_ISIC4       : Sector classification (5 categories)
# YAChiglel17_ISIC4      : Sector classification (17 categories)
# b_year_Group           : Firm age group (5 categories: 0–1, 2–5, 6–10, 11–15, 16+ years)
# ox                     : Ownership type (state-owned, private, foreign-owned, etc.)
# hariutslaga_xelber     : Legal status of the enterprise
# Q2_14_AANT1            : Is the enterprise connected to the e-barimt system? (Yes/No)
# Q2_20_BUGD_AANT1       : Number of employees
# Q4_5_AANT1_q4_4_AANT2_q4_6_AANT3_c1_md1 : Value of assets
# Q4_6_md11_AANT1        : Total revenue
# Q4_7_md1_AANT1         : Cost of intermediate input
# Q4_7_md2_AANT1         : Cost of goods purchased for resale
# Q4_7_md3_AANT1         : Service expense
# Q4_7_md28_AANT1        : Wage expense

# code for creating region (aimag) area df  
mongolia_map <- st_read("data/mng-administrative-divisions-shapefiles/mng_admbnda_adm1_nso_20201019.shp")

df_aimag_code <- data.frame(
  aimag = c(
    13, 14, 15, 12, 16, 17, 18, 19, 21, 20, 22,
    6, 7, 8, 11, 9, 10, 1, 2, 5, 4, 3
  ),
  ADM1_EN = c(
    "Darkhan-Uul", "Dornogovi", "Dundgovi", "Govisumber",
    "Omnogovi", "Selenge", "To'v", "Dornod", "Hentii",
    "Sukhbaatar", "Ulaanbaatar", "Arxangai", "Bayankhongor",
    "Bulgan", "Khovsgol", "Orkhon", "Ovorkhangai",
    "Bayan-Olgii", "Govi-Altai", "Khovd", "Uvs", "Zavkhan"
  ),
  stringsAsFactors = FALSE
)

mongolia_map <- mongolia_map %>%
  left_join(df_aimag_code, by = "ADM1_EN")

area_aimag_df <- mongolia_map %>%
  st_drop_geometry() %>%
  select(aimag, Shape_Area)

area_aimag_df$aimag <- as.character(area_aimag_df$aimag)

##### Loading the main data 

df <- read.spss("data/Data_AANT2021.sav", to.data.frame = TRUE, use.value.labels = F)
UB_dsitrict_code <- read.csv("data/AANT1_UB_Dist_Soumcode.csv")
df <- df %>% left_join(UB_dsitrict_code)
df <- df %>% 
  mutate(
    soum_code = if_else(
      is.na(soum_code),
      NA_character_,
      sprintf("%02d", as.numeric(soum_code))
    ),
    aimag = as.character(aimag),
    code = if_else(
      is.na(soum_code),
      aimag,
      paste0(aimag, soum_code)
    )
  )

# selection of variables form the enterprise census 2021 data
new_df <- df %>%
  filter(Type == "AANT1 ") %>%
  mutate(interm_cost = rowSums(across(c("Q4_7_md1_AANT1", "Q4_7_md2_AANT1", 
                                        "Q4_7_md3_AANT1")), na.rm = TRUE),
         value_added = Q4_6_md11_AANT1 - interm_cost ) %>% 
  filter(interm_cost>0 & value_added >0 & Q2_20_BUGD_AANT1>0 ) %>% 
  mutate(lnY = log(Q4_6_md11_AANT1),
         lnVA = log(value_added),
         lnLC = log(Q4_7_md28_AANT1), # Q2_20
         lnL = log(Q2_20_BUGD_AANT1 ), # Q2_20
         lnK = log(Q4_5_AANT1_q4_4_AANT2_q4_6_AANT3_c1_md1),
         lnM = log(interm_cost))

# Dropped observations below the 1st and above the 99th percentile for each sector and region
vars_to_check <- c("lnY", "lnL", "lnK", "lnM")

df_tfp_clean <- new_df %>%
  filter(lnK >= 0) %>%
  group_by(code, YAChiglel5_ISIC4) %>%
  filter(if_all(all_of(vars_to_check),
                ~ . >= quantile(., 0.01, na.rm = TRUE) &
                  . <= quantile(., 0.99, na.rm = TRUE)))

# =================================================================================
# Estimation TFP using Cobb-Douglas function 
# =================================================================================

fe_model <- feols(lnY ~ lnK + lnL + lnM | YAChiglel5_ISIC4,
                  data = df_tfp_clean)

df_tfp_clean <- df_tfp_clean %>%  ungroup() %>%
  mutate(log_TFP = resid(fe_model))

# =================================================================================
# Creating agglomeration indicators (Localization, Urbanization, Peer TFP)
# =================================================================================
# s for sector, z for region, loc for Localization, urb for Urbanization, avg_prod for peer TFP 
# empl for Employee count, n for Number of firms, _d for density

df_tfp_clean <- df_tfp_clean %>%
  group_by(YAChiglel17_ISIC4,aimag) %>%
  mutate(empl_sz = sum(Q2_20_BUGD_AANT1, na.rm = TRUE),
         n_sz= n())

df_tfp_clean <- df_tfp_clean %>%
  group_by(aimag) %>%
  mutate(empl_z = sum(Q2_20_BUGD_AANT1, na.rm = TRUE),
         n_z= n())

area_aimag_df$aimag <- as.character(area_aimag_df$aimag)
df_tfp_clean <- df_tfp_clean %>% left_join(area_aimag_df , by = "aimag" )

df_tfp_clean <- df_tfp_clean %>%
  mutate(loc = log((empl_sz - Q2_20_BUGD_AANT1)+1),
         urb = log((empl_z-empl_sz)+1),
         loc_d = log((empl_sz - Q2_20_BUGD_AANT1)/Shape_Area+1),
         urb_d = log((empl_z-empl_sz)/Shape_Area+1),
         ln_n_z = log(n_z - n_sz+1),
         ln_n_sz = log(n_sz-1 +1))

df_tfp_clean <- df_tfp_clean %>% group_by(aimag, YAChiglel17_ISIC4) %>% 
  mutate(sum_tfp = sum(log_TFP),
         n= n()) %>% 
  mutate(avg_prod = (sum_tfp-log_TFP)/(n-1)) 

# =================================================================================
# Figure 1: Labor Productivity in Ulaanbaatar
# =================================================================================
# This figure is constructed externally in Excel using the processed data.
# The corresponding file can be found in the "output" folder.
#
# =================================================================================
# Table 1: Average Commute Duration (min) by District
# =================================================================================
# This table is constructed using the Time Use Surveys (TUS) from 2011 and 2023.
# The replication code is provided in a separate script for transparency.
# 
# Note: This step is not required for reproducing the main results of the paper. 
source("scripts/Table1_replication.R")

# =================================================================================
# Figure 2: Spatial distribution of residential areas and urban activity in
# Ulaanbaatar’s core districts
# =================================================================================
# This figure is constructed using Meta’s population data and Overture Maps 
# point-of-interest (POI) data.
#
# The replication code is provided in a separate script for transparency.
#
# Note: This step is not required to reproduce the main results of the paper and 
# may be time-consuming due to data processing requirements. 
source("scripts/Figure2_replication.R")

# =================================================================================
# Table 2: Summary Statistics  
# =================================================================================
# Variables to convert from thousands -> millions
money_vars <- c(
  "Q4_6_md11_AANT1",
  "value_added",
  "Q4_5_AANT1_q4_4_AANT2_q4_6_AANT3_c1_md1",
  "interm_cost"
)

# Variable labels
var_labels <- c(
  Q4_6_md11_AANT1                          = "Gross Revenue (millions of MNT)",
  value_added                              = "Value Added (millions of MNT)",
  Q4_5_AANT1_q4_4_AANT2_q4_6_AANT3_c1_md1 = "Assets (millions of MNT)",
  Q2_20_BUGD_AANT1                         = "# of Employees",
  interm_cost                              = "Intermediate Input Cost (millions of MNT)",
  loc                                      = "Localization",
  urb                                      = "Urbanization",
  avg_prod                                 = "Average Peer Productivity"
)

# Compute descriptive statistics
desc_stats <- df_tfp_clean %>%
  ungroup() %>%
  
  # Temporary conversion only for table creation
  mutate(
    across(all_of(money_vars), ~ .x / 1000)
  ) %>%
  
  select(all_of(names(var_labels))) %>%
  
  pivot_longer(
    cols      = everything(),
    names_to  = "Variable",
    values_to = "Value"
  ) %>%
  
  group_by(Variable) %>%
  
  summarise(
    Mean        = mean(Value, na.rm = TRUE),
    `Std. Dev.` = sd(Value, na.rm = TRUE),
    Median      = median(Value, na.rm = TRUE),
    Min         = min(Value, na.rm = TRUE),
    Max         = max(Value, na.rm = TRUE),
    .groups     = "drop"
  ) %>%
  
  mutate(
    Variable = dplyr::recode(Variable, !!!var_labels)
  )

# Create summary statistics table
desc_stats %>%
  gt() %>%
  
  tab_header(
    title = "Summary Statistics"
  ) %>%
  
  fmt_number(
    columns = c(Mean, `Std. Dev.`, Median, Min, Max),
    decimals = 1
  ) %>%
  
  cols_label(
    Variable    = "Variable",
    Mean        = "Mean",
    `Std. Dev.` = "Std. Dev.",
    Median      = "Median",
    Min         = "Min.",
    Max         = "Max."
  ) %>%
  
  opt_align_table_header("left") %>%
  
  opt_table_font(
    font = "Calibri"
  ) %>%
  
  gtsave("output/table2.pdf")

# =================================================================================
# Figure 3: Localization and Urbanization correlation with TFP estimation
# =================================================================================

# Data preparation
df_plot <- df_tfp_clean %>%
  group_by(aimag, YAChiglel17_ISIC4) %>%
  summarise(
    loc     = mean(loc),
    urb     = mean(urb),
    log_TFP = weighted.mean(log_TFP, value_added),
    .groups = "drop"
  ) %>%
  mutate(
    Region = case_when(
      aimag == 22 ~ "Ulaanbaatar",
      aimag == 13 ~ "Darkhan-Uul",
      aimag == 9  ~ "Orkhon",
      aimag == 16 ~ "Umnugovi",
      TRUE        ~ "Other"
    )
  )
# Localization
p_loc <- ggplot(df_plot, aes(x = loc, y = log_TFP, shape = Region)) +
  geom_point(size = 2) +
  geom_smooth(
    aes(group = 1),         
    method  = "lm",
    formula = y ~ x ,
    se      = TRUE,
    color   = "black",
    size    = 0.7
  ) +
  labs(x = "Localization", y = "Log TFP", title = "Localization vs. TFP") +
  scale_shape_discrete(name = "Region:") +
  theme_minimal() +
  theme(legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))

#  Urbanization
p_urb <- ggplot(df_plot, aes(x = urb, y = log_TFP, shape = Region)) +
  geom_point(size = 2) +
  geom_smooth(
    aes(group = 1),          
    method  = "lm",
    formula = y ~ x ,
    se      = TRUE,
    color   = "black",
    size    = 0.7
  ) +
  labs(x = "Urbanization", y = "Log TFP", title = "Urbanization vs. TFP") +
  scale_shape_discrete(name = "Region:") +
  theme_minimal() +
  theme(legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))

# Combine side-by-side with shared legend at bottom
combined <- (p_loc + p_urb) +
  plot_layout(ncol = 2, guides = "collect") &
  theme(legend.position = "bottom")

# Save as PNG
ggsave(
  filename = "output/Figure3.png",
  plot = combined,
  width = 10,     
  height = 5,
  dpi = 300        
)

# =================================================================================
# Table 3: Estimation of Agglomeration Impacts on Firm Productivity
# =================================================================================

# Fixed effects 
f_effect <- c("YAChiglel17_ISIC4", "b_year_Group", "ox",
              "hariutslaga_xelber", "Q8_1", "Q2_14_AANT1")

# Cluster
cl <- "YAChiglel17_ISIC4"

# Models
ls_01 <- feols(log_TFP ~ avg_prod,
               data = df_tfp_clean,
               fixef = f_effect,
               cluster = cl)

ls_02 <- feols(log_TFP ~ urb,
               data = df_tfp_clean,
               fixef = f_effect,
               cluster = cl)

ls_03 <- feols(log_TFP ~ loc,
               data = df_tfp_clean,
               fixef = f_effect,
               cluster = cl)

ls_04 <- feols(log_TFP ~ avg_prod + urb + loc,
               data = df_tfp_clean,
               fixef = f_effect,
               cluster = cl)

ls_05 <- feols(log_TFP ~ avg_prod + urb + I(urb^2) + loc + I(loc^2),
               data = df_tfp_clean,
               fixef = f_effect,
               cluster = cl)

ls_06 <- feols(log_TFP ~ avg_prod + urb + loc + I(loc^2),
               data = df_tfp_clean,
               fixef = f_effect,
               cluster = cl)

# Report table
cm = c(
  'avg_prod'    = 'Peer TFP',
  'urb'    = 'Urbanization',
  'I(I(urb^2))'    = ' Quadratic term of Urbanization',
  'loc'    = 'Localization',
  'I(I(loc^2))'    = ' Quadratic term of Localization')

note2table1 <- c(
  "Notes: Each regression includes fixed effects for 17 sectors, firm age groups, ownership types, legal forms, R&D investment, and technology adoption (E-barimt). Standard errors, clustered at the sector level, are reported in brackets. *** Significant at the 1% level, ** Significant at the 5% level, * Significant at the 10% level.")
titletable1 <- "Table 3: Estimation of Agglomeration Impacts on Firm Productivity"

msummary(
  list(ls_01, ls_02, ls_03, ls_04, ls_05, ls_06),
  stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  estimate = "{estimate}",
  statistic = "[{std.error}]{stars}",
  coef_map = cm,
  gof_omit = 'AIC|BIC|RMSE|Within|Std.Errors|FE',
  notes = note2table1,
  title = titletable1,
  output = "output/table3.txt"
  )

# =================================================================================
# Figure 4: Inverted U-shaped Curve of Over-Agglomeration
# =================================================================================

# prepare the data aimag level averages
table1 <- df_tfp_clean %>% group_by(aimag) %>% 
  summarise(n = n(), 
            mean = weighted.mean(log_TFP, value_added),
            urb = weighted.mean(urb, value_added),
            loc =weighted.mean(loc, value_added),
            n_z = sum(n_z),
            n_sz = sum(n_sz))


# Define the quadratic function for the model (6)
coef_lin  <- ls_06$coefficients[3]
coef_quad <- ls_06$coefficients[4]
f <- function(x) coef_lin * x + coef_quad * x^2

x_vals <- seq(0, 15, by = 0.1)
curve_df <- tibble(
  loc = x_vals,
  tfp = f(x_vals))

x_max <- -coef_lin / (2 * coef_quad)
y_max <- f(x_max)

# Adding aimag names
table1 <- table1 %>%
  mutate(
    Region = case_when(
      aimag == 22 ~ "Ulaanbaatar",
      aimag == 13 ~ "Darkhan-Uul",
      aimag == 9  ~ "Orkhon",
      aimag == 16 ~ "Umnugovi",
      aimag == 17 ~ "Selenge",
      TRUE        ~ "Other"
    )
  )

# Plot with ggplot2
fig4 <- ggplot() +
  geom_line(data = curve_df, aes(x = loc, y = tfp),
            color = "grey40", size = 1) +
  geom_point(data = table1,
             aes(x = loc, y = f(loc), shape = Region),
             size = 3, alpha = 0.9) +
  geom_vline(xintercept = x_max,
             linetype = "dashed", color = "grey40") +
  labs(
    title    = "Agglomeration Effect on TFP",
    subtitle = "Fitted quadratic trend with regional observations",
    x        = "ln(Localization)",
    y        = "ln(TFP)",
    color    = "Region"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    plot.title      = element_text(face = "bold"),
    plot.subtitle   = element_text(size = 12)
  ) +
  scale_color_brewer(palette = "Set1") 
 
#saving results as a picture
ggsave(
    filename = "output/Figure4.png",
    plot = fig4,
    width = 10,
    height = 5,
    dpi = 300,
    bg = "white"
  )

# ================================================================================================
# Table 4: Fixed-Effects Estimation of Agglomeration Impacts on Firm Productivity by Sector
# ================================================================================================

# Model estimation
# In the ISIC 5-sector classification, agriculture and construction represents only one sector 
# within the ISIC 17-sector classification. Therefore, to obtain clustered 
# standard errors for other 3 sectors, the first nd third regression is estimated separately from the loop below.

temp_df <- df_tfp_clean %>% filter(YAChiglel5_ISIC4 == 1)
ls_11 <- feols(log_TFP ~ avg_prod + urb + loc + I(loc^2),
               data = temp_df,
               fixef = f_effect)


temp_df <- df_tfp_clean %>% filter(YAChiglel5_ISIC4 == 3)
ls_13 <- feols(log_TFP ~ avg_prod + urb + loc + I(loc^2),
               data = temp_df,
               fixef = f_effect)


for (i in c(2, 4:5)) {
  # Filter data
  temp_df <- df_tfp_clean %>% filter(YAChiglel5_ISIC4 == i)
  
  # Fit model
  model <- feols(
    log_TFP ~ avg_prod + urb + loc + I(loc^2),
    data   = temp_df,
    fixef  = f_effect,
    cluster = cl
  )
  assign(paste0("ls_1", i), model)
}

# Report table
titletable4 <- "Table 4: Fixed-Effects Estimation of Agglomeration Impacts on Firm Productivity by Sector"
note2table4 <- c("(1) Agriculture, (2) Manufacturing, (3) Construction, (4) Trade, (5) Service")

msummary(
  list(ls_11, ls_12, ls_13, ls_14, ls_15),
  stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
  estimate = "{estimate}",
  statistic = "[{std.error}]{stars}",
  coef_map = cm,
  gof_omit = 'AIC|BIC|RMSE|Within|Std.Errors|FE',
  notes = note2table4,
  title = titletable4,
  output = "output/table4.txt"
)
 
 # ================================================================================================
 # Table A1: Fixed-Effects Estimation of Agglomeration Impacts on Firm Productivity by Firm Size
 # ================================================================================================

 # Model estimation
 for (i in 1:3) {
   # Filter data
   temp_df <- df_tfp_clean %>% filter(Emp_Total_AANT1_AANT2_AANT3_Group == i)
   
   # Fit model
   model <- feols(
     log_TFP ~ avg_prod + urb + loc + I(loc^2),
     data   = temp_df,
     fixef  = f_effect,
     cluster = c("YAChiglel17_ISIC4")
   )
   assign(paste0("ls_2", i), model)
 }
 
 temp_df <- df_tfp_clean %>% filter(Emp_Total_AANT1_AANT2_AANT3_Group == 4)
 ls_24 <- feols(
   log_TFP ~ avg_prod + urb + loc,
   data   = temp_df,
   fixef  = f_effect,
   cluster = c("YAChiglel17_ISIC4")
 )
 
 # Report table
 titletableA1 <- "Table 2: Fixed-Effects Estimation of Agglomeration Impacts on Firm Productivity by firm size"
 note2tableA1 <- c("(1) Micro (1-9), (2) Small (10-19), (3) Medium (20-49), (4) Large (50+)")
 
 
 msummary(
   list(ls_21, ls_22, ls_23, ls_24),
   stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
   estimate = "{estimate}",
   statistic = "[{std.error}]{stars}",
   coef_map = cm,
   gof_omit = 'AIC|BIC|RMSE|Within|Std.Errors|FE',
   notes = note2tableA1,
   title = titletableA1,
   output = "output/tableA1.txt"
 )
 
 # ================================================================================================
 # Table A2: Robustness Check: Alternative Agglomeration Measures (Firm Counts)
 # ================================================================================================
 
 # Model estimation
 ls_31 <- feols(log_TFP ~ avg_prod,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_32 <- feols(log_TFP ~ ln_n_z,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_33 <- feols(log_TFP ~ ln_n_sz,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_34 <- feols(log_TFP ~ avg_prod + ln_n_z + ln_n_sz,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_35 <- feols(log_TFP ~ avg_prod + ln_n_z + I(ln_n_z^2) + ln_n_sz + I(ln_n_sz^2),
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_36 <- feols(log_TFP ~ avg_prod + ln_n_z + ln_n_sz + I(ln_n_sz^2),
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 # Report table
 titletableA2 <- "Table A2: Robustness check: Alternative Agglomeration Measures (Firm Counts)"
 
 cm = c(
   'avg_prod'    = 'Peer TFP',
   'ln_n_z'    = 'Urbanization',
   'I(I(ln_n_z^2))'    = ' Quadratic term of Urbanization',
   'ln_n_sz'    = 'Localization',
   'I(I(ln_n_sz^2))'    = ' Quadratic term of Localization')
 
 msummary(
   list(ls_31, ls_32, ls_33, ls_34, ls_35,  ls_36),
   stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
   estimate = "{estimate}",
   statistic = "[{std.error}]{stars}",
   coef_map = cm,
   gof_omit = 'AIC|BIC|RMSE|Within|Std.Errors|FE',
   notes = note2table1,
   title = titletableA2,
   output = "output/tableA2.txt"
 )
 
 # ================================================================================================
 # Table A3: Robustness Check: Dependent Variable as Labor Productivity
 # ================================================================================================
 
# Prepare the data 
 df_tfp_clean <- df_tfp_clean %>% mutate(lnLP = log(value_added / Q2_20_BUGD_AANT1))
 df_tfp_clean <- df_tfp_clean %>% group_by(aimag, YAChiglel17_ISIC4) %>% 
   mutate(sum_lp = sum(lnLP),
          n= n()) %>% 
   mutate(avg_lp = (sum_lp-lnLP)/(n-1))
 
 # Model estimation
 ls_41 <- feols(lnLP ~ avg_lp,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_42 <- feols(lnLP ~ urb,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_43 <- feols(lnLP ~ loc,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_44 <- feols(lnLP ~ avg_lp + urb + loc,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_45 <- feols(lnLP ~ avg_lp + urb + I(urb^2) + loc + I(loc^2),
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_46 <- feols(lnLP ~ avg_lp + urb + loc + I(loc^2),
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 # Report table
 titletableA3 <- "Table A3: Robustness check -  Dependent Variable as Labor Productivity"
 
 cm = c(
   'avg_lp'    = 'Average LP',
   'urb'    = 'Urbanization',
   'I(I(urb^2))'    = ' Quadratic term of Urbanization',
   'loc'    = 'Localization',
   'I(I(loc^2))'    = ' Quadratic term of Localization')
 
 msummary(
   list(ls_41, ls_42, ls_43, ls_44, ls_45,ls_46),
   stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
   estimate = "{estimate}",
   statistic = "[{std.error}]{stars}",
   coef_map = cm,
   gof_omit = 'AIC|BIC|RMSE|Within|Std.Errors|FE',
   notes = note2table1,
   title = titletableA3,
   output = "output/tableA3.txt"
 )
 
 # ================================================================================================
 # Table A4: Robustness Check: Urbanization and Localization as Distance-Based Metrics
 # ================================================================================================
 
 # Model estimation
 ls_51 <- feols(log_TFP ~ avg_prod,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_52 <- feols(log_TFP ~ urb_d,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_53 <- feols(log_TFP ~ loc_d,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_54 <- feols(log_TFP ~ avg_prod + urb_d + loc_d,
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_55 <- feols(log_TFP ~ avg_prod + urb_d + I(urb_d^2) + loc_d + I(loc_d^2),
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 ls_56 <- feols(log_TFP ~ avg_prod + urb_d + loc_d + I(loc_d^2),
                data = df_tfp_clean,
                fixef = f_effect,
                cluster = cl)
 
 # Report table 
 cm = c(
   'avg_prod'    = 'Peer TFP',
   'urb_d'    = 'Urbanization',
   'I(I(urb_d^2))'    = ' Quadratic term of Urbanization',
   'loc_d'    = 'Localization',
   'I(I(loc_d^2))'    = ' Quadratic term of Localization')
 
 titletableA4 <- "Table A4: Robustness check: Urbanization and Localization as Distance-based Metrics"
 
 msummary(
   list(ls_51, ls_52, ls_53, ls_54, ls_55, ls_56),
   stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01),
   estimate = "{estimate}",
   statistic = "[{std.error}]{stars}",
   coef_map = cm,
   gof_omit = 'AIC|BIC|RMSE|Within|Std.Errors|FE',
   notes = note2table1,
   title = titletableA4,
   output = "output/tableA4.txt"
 )